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ABSTRACT 


This  report  is  a  brief  survey  of  some  currently  used  techniques  of  periodic 
analysis  of  sampled  data.  Almost  every  periodic  analysis  of  empirical  data  eventually 
relies  on  one  of  the  following  three  techniques: 

1  .  Regression  analysis  on  a  trigonometric  polynomial 

2.  Periodogram  analysis 

3.  The  estimation  of  power  spectra. 

Each  of  these  techniques  is  best  applied  under  distinct  modelling  assumptions,  and  an 
attempt  has  been  made  to  discuss  their  applicability  and  their  range  of  validity  under 
various  circumstances. 

After  a  brief  introduction  in  Chapter  I,  Chapter  II  proceeds  to  describe  the  models 
in  use  and  the  statistical  properties  of  the  estimates  associated  with  them  The  effects 
of  sampling  and  some  of  their  implications  in  practical  situations  are  treated  in  Chapter  III. 
Chapter  IV  relates  several  common  question  areas  to  the  theory  developed  in  the  earlier 
portions  of  the  report,  and  presents  some  illustrations  of  the  power  of  the  techniques 
discussed.  Finally,  Chapter  V  gives  the  briefest  outline  of  the  general  requirements  of 
a  comprehensive  program  in  time  series  data  analysis.  Numbers  in  square  brackets 
refer  to  entries,  in  alphabetic  order,  in  the  appended  bibliography. 
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I. 


INTRODUCTION 


In  the  processing  of  empirical  data,  a  frequency  analysis  is  often  required  as  a 
partial  characterization  of  the  process  generating  the  data.  In  particular,  one  hopes, 
on  the  basis  of  such  an  analysis,  to  draw  some  inference  on  the  presence  of  periodicities 
in  the  underlying  process.  Two  approaches  to  frequency  analysis  (harmonic  analysis  and 
spectral  estimation)  are  commonly  taken,  sometimes  without  full  awareness  of  their 
implication  and  range  of  validity.  The  sources  of  ambiguity  are  several: 

1.  There  may  be  a  basic  ignorance  about  the  nature  of  the  underlying 
phenomenon,  in  which  case  the  choice  of  a  suitable  model  becomes  the 
main  stumbling  block.  Indeed,  the  following  kinds  of  alternatives  have 
to  be  considered,  among  others: 

-  the  phenomenon  may  be  periodic,  in  which  case  all  frequencies  are 
multiples  of  a  single  base  frequency; 

-  the  phenomenon  may  be  aperiodic,  in  which  case  frequencies  of 
interest  may  be  either  discrete  or  continuous ; 

-  the  phenomenon  may  be  thought  of  as  a  single  unique  occurrence, 
or  it  may  be  taken  as  one  of  many  possible  similar  realizations. 

2.  Each  of  these  alternatives  represents  a  different  point  of  view  and  requires 
different  techniques  in  the  analysis  of  data  and  the  interpretation  of  results. 

In  particular,  there  is  an  apparent  formal  similarity  between  conventional 
techniques  of  Fourier  Analysis  and  the  estimation  of  power  spectra,  which 
leads  to  the  temptation  of  using  one  for  the  other  although  they  apply  to 
distinctly  different  models  and  have  essentially  different  statistical  prop¬ 
erties.  This  gives  rise  to  modelling  uncertainty  and  the  attendant  problems 
in  interpretation  of  results. 

3.  The  current  literature  on  the  subject  of  frequency  analysis  of  empirical 
time  series  draws  its  techniques  and  terminology  from  at  least  three  fields: 

-Classical  Statistical  Theory,  [3,  11,  12,  16] 

-Purely  mathematical  disciplines  such  as  approximation  theory,  [1] 
representation  theory,  [6]  Fourier  Analysis  and  related  fields, 

[19,  25,  29] 
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-The  frequency  analysis  of  linear  systems  in  communications 
engineering,  [14,  19,  29]. 

Here,  unfamiliarity  with  one  or  another  language  is  the  stumbling  block. 

This  report  has  been  written  in  an  attempt  to  elucidate  the  significant  relations 
governing  the  applicability  and  validity  of  the  theory.  The  subject  to  be  covered  is  large 
and  growing,  so  that  our  presentation  is  anything  but  complete.  Moreover,  the  theory 
requires  some  fairly  advanced  mathematical  tools;  therefore  the  presentation  is  primarily 
discursive.  Every  argument  can  however  be  rigorously  demonstrated  without  alteration 
of  content. 

We  shall  be  considering  primarily  zero  mean  second  order  stationary  stochastic 
processes,  as  there  exists  hardly  any  theory  for  evolutive*  time  series.  Almost  all 
references  dealing  with  the  latter  generally  convert  them  to  something  resembling  station - 
arity  by  means  of  some  elementary  statistical  device. 

Chapter  II  describes  the  models  which  are  most  commonly  employed  in  periodic 
analysis.  These  are  a  cyclic  periodic  model,  a  cyclic  aperiodic  model  and  an  oscillatory 
model.  The  cyclic  models  are  presented  primarily  because  they  offer  a  good  background 
for  the  treatment  of  the  more  sophisticated  oscillatory  model.  Emphasis  will  be  placed 
on  the  statistical  properties  of  the  various  estimation  techniques  presented. 

Chapter  III  will  treat  several  consequences  of  sampling  and  the  implications  which 
follow  in  practical  applications.  In  particular  the  phenomena  of  folding  and  aliasing  will 
be  described. 

Chapter  IV  has  a  two  fold  purpose;  on  the  one  hand  several  common  question  areas 
will  be  elucidated  in  the  light  of  the  foregoing  theory;  on  the  other,  several  illustrations 
will  be  given  of  the  power  of  the  techniques  described  in  the  earlier  chapters. 

Finally,  Chapter  V  outlines  briefly  the  necessary  data  processing  and  analysis 
areas  which  should  be  included  in  a  comprehensive  program  of  periodic  analysis  of 
empirical  data. 

*An  evolutive  time  series  is  one  with  a  non  zero  time  varying  mean. 
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II.  MODELS 


We  shall  examine  three  models  which  find  most  frequent  application.  The  first  two 
are  cyclic  models;  they  are  the  simpler  of  the  three  and  have  been  studied  extensively  in 
the  last  forty  years.  They  are  based  on  the  assumption  that  the  time  series  to  analyze  is 
the  sum  of  a  deterministic  component  and  white  noise.  The  deterministic  component  is  a 
finite  sum  of  sine  and  cosine  terms  whose  frequencies  may  or  may  not  be  multiples  of  a 
single  base  frequency.  The  third  model,  known  as  oscillatory,  is  based  on  the  assumption 
that  the  time  series  is  a  realization  of  a  stochastic  process  possessing  a  continuous  power 
spectrum.  A  heuristic  difference  between  cyclic  and  oscillatory  models  is  that  for  cyclic 
processes,  different  realizations  are  expected  to  look  alike  except  for  small  fluctuations, 
whereas  two  single  realizations  of  an  oscillatory  process  may  look  entirely  different  from 
each  other.  Another  way  of  stating  this  is  to  say  that  different  realizations  of  a  cyclic 
process  agree  in  mean  whereas  different  realizations  of  an  oscillatory  process  agree  in 
quadratic  mean.  While  the  first  two  models  rest  mainly  on  the  application  of  classical 
techniques,  the  third  model  has  only  been  developed  vigorously  during  the  last  fifteen  years 
and  is  still  the  subject  of  current  research. 

A  fourth  model  is  of  great  practical  interest  and  deals  with  processes  which  are  a 
mixture  of  cyclic  and  oscillatory  processes;  it  assumes  a  superposition  of  oscillatory  and 
cyclic  processes  and  a  spectrum  which  need  not  be  continuous.  Analysis  of  this  kind  of 
model  is  almost  non-existent  in  the  literature  but  it  appears*  that  work  is  currently  in 
progress  on  spectral  estimates  which  converge  to  the  derivative  of  the  integrated  spectrum 
where  it  exists,  and  to  the  magnitude  of  the  jump  where  continuity  fails.  This  kind  of  model 
shall  not  be  considered  in  this  report. 

In  general,  a  realization  of  a  continuous  process  {x^}  will  be  represented  as  x(t).  A 
discrete  process  {x  }  will  have  a  realization  x^.  In  the  present  chapter,  we  shall  not  distin¬ 
guish  between  a  realization  of  a  discrete  process  and  a  sample  of  a  realization  of  a  continuous 
process;  both  will  be  written  as  x^.  For  the  present  this  is  a  minor  distinction  although  we 
shall  see  later  precisely  under  what  conditions  it  becomes  important. 

A.  First  and  Second  Model:  Cyclic  Structure 

We  shall  assume  in  this  case  that  the  time  series  x^  is  generated  by  a  process 

of  the  form: 

*Doctoral  Dissertion  in  progress  under  E.  Parzen.  (Private  communication  by  L.  Gardner). 
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X  = 
t 


■I 


(  a.  cos  A.t  +  b.  sin  X.t  )  +  £ 

J  J  J  J  t 


j=l 


where*  0  <  rr ,  a.  and  b.  are  real  and  £  is  Gaussian,  uncorrelated,  with  mean  zero 


J 


J 


J 


We  can  assume  that  the  frequencies  A  are  all  multiples  of  a  base 


and  variance  cr  .  There  are  two  possibilities  open  at  this  point: 

(1) 

frequency  known  a  priori. 

(2)  We  know  neither  the  frequencies  nor  the  amplitudes  of  the  process 
and  are  trying  to  estimate  both,  and  in  particular,  the  frequency 
associated  with  the  peak  amplitude. 

The  first  possibility  leads  to  a  straightforward  application  of  linear  regression  theory 
[  l6,Vol  II,  pp.  141-174]  with  trigonometric  functions  substituted  for  polynomials,  and 
is  equivalent  to  a  least  squares  fit  of  an  empirical  function  by  a  trigonometric  polynomial 
[  34,  pp.  14-19]  .  The  second  possibility  leads  to  the  calculations  of  the  periodogram  and 
is  commonly  referred  to  as  "Harmonic  Analysis"  of  empirical  data. 

1.  First  Interpretation:  Frequencies  Known,  Amplitudes  Unknown 
[  25,  Ch  I] 

In  this  case,  the  process  is  a  superposition  of  white  noise  and  a 
periodic  deterministic  component.  All  frequencies  are  multiples  of  a  base  frequency 
and  we  can  write 


t  =  a  +  )  (a 

t  0  U  v. 

j=l  J 


2nv. 


2rrv, 


cos 


t  +  /?  sin 

v 

j 


t  >  +  «,. 


*  A.  is  restricted  to  be  less  than  tt  because  in  the  concrete  situation  of  a  sample  of 
J  =  .  .  .  ...  .  ,  T 


size  T,  we  have  X 


T  and  j  cannot  be  allowed  to  exceed  -  ,  (resulting  in  an 


upper  bound  of  tt  for  A.)  lest  there  be  more  than  T  coefficients  to  evaluate  on  the  basis 
of  a  sample  of  size  T. 
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We  assume  m  =  x  -  |  to  be  periodic  with  period  p,  m  ,  =  m  ,  and  further,  that 

t  t  t  K  r  r,  t+kp  t> 

the  number  of  observations  T  is  a  multiple  (large  usually  )  of  the  period  p.  |  is  white 
noise  as  previously  described.  The  v  ,  j=l— >q  are  a  subset  of  the  sequence  1, 2. . .  ,  p 
and  are  used  to  indicate  that  not  every  harmonic  need  be  represented*.  What  is  wanted 


is  an  estimate  of  the  coefficients  a 


v  ,  Bv  .  Since  £  is  Gaussian,  the  least  squares 
j  J  1 


estimate  is  the  maximum  likelihood  estimate  of  a  and  B  and  we  find  that 


t=l 


V  =  T 

J 


2 TTV , 

X  cos  - J 

t  p 


t=l 


B 


=  4  I 


2  TTV. 


x  sin 
t 


t=l 


t  ,  j  =  1— *q 


are  the  estimates  of  the  d s  and  (3 's. 

Furthermore,  the  A  and  B  are  independent  normal  variables  with 

"j  2a2 

means  a  and  B  p  respectively,  and  their  variances  are  .which  go  to  0  as  sample 

V)  j 

size  T— ►°°. 


Setting  the  squared  amplitude  at  a  given  frequency  to 

2  2  2 

R  =  A  +  B  , 

V  V  V 

j  j  J 

2  2  2  2 
it  can  be  shown  that  R  is  an  estimate  of  the  true  squared  amplitude  p(  ~  av  +  P 

V)  j  j  j 

with  the  properties 


*  We  have  chosen  the  v  scheme  of  indexing  the  a's  and  B  's  rather  than  letting  j  =  1 — ►p 
and  allowing  some  a^'s  and  b/s  to  be  zero  because  the  presence  of  the  noise  ^  will  insure 
that  none  of  the  estimated  oi's  and  B' s  will  be  zero,  whereas  we  might  know  a  priori  that 
some  of  them  are. 
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2 

P  + 
v. 

J 


„  _  2  8a' 

Var  R  =  T 

J 


2  r 


2  2a 

V+  *T 

J 


In  other  words,  R  is  a  consistent  asymptotically  unbiased  estimate*  of  the  squared 

j 

amplitude  of  the  v  th  harmonic. 

It  may  be  pointed  out  that  very  rarely  is  there  enough  information  to 

assume,  in  applications,  that  the  fundamental  frequency  is  known.  The  analysis  of  this 

model  is  thus  primarily  a  theoretical  exercise  in  regression  analysis. 

2.  Second  Interpretation:  Amplitude  and  Frequencies  Unknown 

This  model,  as  the  previous  one,  is  a  finite  parameter  model  of  the 

form  q 

x  =  /  (  q;.  cos  o>.t  +  /3 .  sin  oj.ti+f  ,  t  =  1— -  T  . 
t  Lj  \  J  J  J  J  /  t 

j=l 


The  constant  term  has  been  omitted  as  it  can  be  estimated  and  removed.  Again,  £  is 
white  Gaussian  noise;  the  only  restriction**  on  the  o>  is 

j 

0  <  ay.  <  7r  ; 


no  longer  are  the  uk  harmonically  related.  The  problem  this  time  is  to  extract  the 
periods  27r/uK  from  the  sequence  x^. 

The  technique  of  the  periodogram  proceeds  as  follows : 


*Unbiased  A  the  mean  value  of  the  estimate  equals  the  true  value. 

Consistent  A  the  estimate  converges  in  probability  to  the  true  value.  Sufficient  for 
consistency  is  that  Var(estimate)  -*0  as  sample  size  -*  °°. 

Asymptotically  A  . . .  as  the  sample  size  goes  to  °° . 

The  symbol  A  is  to  be  read  as:  "means  by  definition". 

**cf.  footnote  page  4, 
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Calculate  the  periodogram 


at  all  points 


~  2ttT 


V  icjt 

)  x  -e 


t  =  1 


[A2(o;)+B2(a;)] 


2nV  ,  T  . 

=  -Tjr-  ,  f=  1"*2_1’  (evenT)> 


where 


A(w)=-  N  cos  wt 
t=  1 
T 

B(w)  =  |^  xt  sin  cot 
t  =  1 


are  the  Fourier  coefficients  of  the  data  relative  to  a  trial  frequency  co.  In  short,  x  is 
subjected  to  an  ordinary  Fourier  Analysis  as  if  the  noise  term  £  were  absent. 

Now,  it  is  a  characteristic  property  of  the  periodogram  <!»  (co)  that  for 
large  samples  it  is  large  whenever  co  is  equal  to  one  of  the  frequencies  or  and  that  it  is 
small  whenever  co  f  or  [11,  p.  92]*.  Thus,  as  one  expects  large  peaks  (for  large  samples) 
in  4>^,(co)  whenever  co  =  co.,  it  is  natural  to  say: 

"Let  us  calculate  ^(co)  at  all  frequencies  of  interest  and  single  out 
the  peaks  in  ^(co),  inferring  therefrom  that  a  frequency  component 
is  predominantly  present  in  the  process  x^_  wherever  a  peak  occurs 
in  $T(co).  " 

The  question  of  the  statistical  significance  of  the  peaks  in  the  periodo¬ 
gram  arises  immediately.  Several  tests  of  significance  have  been  developed  and  are 
briefly  described  in  [15] .  Ideally,  one  should  like  to  state  that  peaks  occur  if  and  only 


th 


*Heuristically,  the  periodogram  is  the  square  of  the  amplitude  c  of  the  v  component; 

if  we  write  the  v^1  component  as  c^  cos  (ctjt  +  <p  ).  It  is  known  that  when  £  =  0  , 

T 

1  V  2  1 

T_Xt  2 
t=  1  v 

so  that  c^  is  proportional  to  the  contribution  of  the  v  component  to  the  total  energy 
dissipated  by  the  process  (cf.  Parseval's  relation  and  completeness). 
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if  a  corresponding  frequency  is  present  in  the  signal;  this  however  is  not  so,  as  can 
be  directly  inferred  from  the  degradation  imposed  by  the  noise  term  We  shall 
mention  only  one  test,  due  to  Fisher,  which  is  independent  of  sample  size. 

To  assess  the  significance  of  the  largest  peak  in  the  periodogram 
compute  the  statistic* 


g=  <f>  (o>  ) 

T  max 


(T+l) 


$T(ax) 


j=1 


where  co  is  the  frequency  corresponding  to  the  maximum  amplitude  of  $  (cu),  where 
max  ~  .  i 

27 r  i 

T  is  odd  and  ca.  =  •  The  distribution  of  g  can  be  explicitly  calculated.  Details  of 

usage  are  given  in  Ref.  [11,  pp.  91  ff] .  Other  tests  are  referred  to  in  Ref.  [15] . 


The  distribution  of  the  values  of  the  periodogram  of  a  set  of  T  random, 

2  o*2  2 

uncorrelated,  normal  (0,  a  )  numbers  is-^-  X ^  e*Pectecl  value  of  the  periodogram 

<j2  <j^ 

is  —  ,  its  variance  is— — j  .  The  probability  that  the  periodogram  exceed  k  times  its 


47 r ' 


expected  value  is 


p^T(“>>k27 


} 


} 


-k 

=  e 


This  is  the  basis  of  Shuster's  original  idea  of  using  the  periodogram  to  reject  the  hypo¬ 
thesis  -  whenever  $^(0;)  is  large  -  that  no  periodic  component  is  present  in  the  data  at 
frequency  An  extensive  exposition  of  periodogram  analysis  is  given  by  Bartels  in 
Ref.  [  15  ]  and  Stumpf  in  Ref  [25] 
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B.  Third  Model:  Oscillatory  Processes 

In  very  few  cases  is  it  reasonable  on  a  priori  grounds  to  decide  that  a  process 
is  of  the  finite  parameter  type.  The  previous  section  briefly  described  two  such  models: 

-a  periodic  model  where  all  frequencies  were  an  integral  multiple  of  a  base 
frequency,  sampled  over  an  interval  containing  several  periods; 

-an  aperiodic  model  where,  while  discrete,  the  spectral  components  were 
no  longer  harmonically  related. 

We  will  now  consider  the  following: 

-an  aperiodic  model  with  a  continuous  spectrum. 

The  cornerstone  of  the  development  is  the  well  known  theorem  which  asserts  that  for  a 
wide  class  of  stochastic  processes  the  autocorrelation  function  and  the  power  spectral 
density  are  completely  equivalent  descriptions  of  the  second  order  properties  of  the 
process. 

a 

1.  Preliminaries :  In  general,  if  f(co  )  is  an  estimate  of  f(to  ),  we  shall 

a 

use,  as  an  evaluation  criterion  of  f,  the  mean  square  error 

r\ 

MSEf=  Bias  f+Varf 

where 

Bias  f  =  (5  (f  -  f ) 

Var  f  =  (f  (f  -(ff)2 

There  are  several  other  criteria  developed  in  the  literature  [15,  18,  31] .  We  shall  not 
pursue  these. 

We  recall  that 

a. 

f  is  unbiased  if  Bias  f  =  0 

A  A 

f  is  asymptotically  unbiased  if  Bias  f  — -  0  as  sample  size  T  — -  «> 

A  A. 

f  is  called  consistent  if  it  is  unbiased  and  Var(f)  — ■  0  as  T  — - 
Sometimes,  in  the  sequel,  we  shall  write  f(ajQ)  instead  of  f(u>).  This  is  done  to  stress 
the  fact  that  we  are  concerned  with  pointwise  estimates,  i.  e. ,  we  are  estimating  a 
single  value,  the  value  that  f  takes  at  w  .  We  are  not  evaluating  a  function.  The  esti¬ 
mation  of  the  function  is  done  pointwise,  one  by  one,  in  practical  estimations.  In  this 
chapter  the  emphasis  is  on  statistical  properties  of  estimates. 
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We  shall  be  concerned  with  zero  mean  second  order  stationary  processes 
having  continuous  spectra  and  possessing  no  periodic  components  corresponding  to  spectral 
lines*.  Roughly  speaking,  a  stochastic  process  is  a  collection  of  functions  {x^}.  Each 
function  is  a  possible  realization  of  the  process  and  it's  occurrence  is  thought  to  be 
governed  by  a  probability  law.  A  zero  mean  process  obeys  the  relation: 

fx  =  0  for  all  t. 
t 

Stationarity  implies  that  the  probabilities  governing  the  occurrence  of  certain  values 


for  any  n,  are  independent  of  the  time  origin,  and  that  these  values  are  thus  just  as 
likely  to  occur  at  one  or  the  other  of  the  two  time  sequences 


V  t2’ 


t 

n 


t  .+  h,  t^+h,  .  .  . ,  t  +h, 

12  n 

regardless  of  the  shift  h.  Second  order  stationarity  means  that  all  moments  up  to 
moments  of  second  order  are  stationary.  In  particular, 


<f 


X  X 


;t+  Ik  I 


=  r(k) 


depends  only  on  the  lag  k.  The  function  y(k)  is  the  autocovariance  function  of  the  process. 

2 

It  is  sometimes  normalized  by  dividing  by  y(0)  =  a  and  the  resulting  function  is  the  auto¬ 
correlation  function.  A  Gaussian  process  is  completely  specified  by  its  first  and  second 
order  properties.  Most  of  the  theory  is  developed  for  Gaussian  processes.  If  the  Gaussian 
assumption  is  dropped,  results  pertaining  to  expected  values  are  not  affected;  but  results 
pertaining  to  the  variability  of  estimates  become  approximations. 


*Processes  with  mixed  spectra  have  been  very  scantily  treated.  One  reference  is  [30] ; 
also,  see  footnote  on  page  3. 
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There  are  two  pairs  of  alternative  ways  of  presenting  results,  corres 
ponding  to  the  real-complex  presentation  on  the  one  hand,  and  the  continuous -discrete 
cases  on  the  other.  The  first  dichotomy  is  purely  a  matter  of  convenience  in  notation. 
The  second  corresponds  to  the  distinction  between  Fourier  series  and  Integrals  (the 
Fourier  series'  coefficients  are  the  Fourier  transform  of  a  periodic  function).  We 
shall  stay  with  the  discrete  case  as  this  is  the  one  most  frequently  encountered  in 
practice. 


The  basic  theorem  states  that  under  rather  general  conditions  insuring 

that  the  covariance  function  y(k)  =  £  x  x  =  y(-k)  goes  to  zero  sufficiently  rapidly  as 

t  t  i  K 

the  lag  goes  to  infinity,  the  second  order  zero  mean  stationary  process  {x^}  possesses 
a  spectrum 


oo 


k=-oo 


which  is  the  Fourier  transform  of  the  autocovariance  y(k),  and  the  converse  relation, 

7 r 

y(k)  =  jVkaJf(o;)du;,  k  =  0,  +  l,+2,.  .  . 

-7T 

also  holds.  (This  is  the  discrete  version  of  the  Wiener -Khintchin  theorem  [34,  pp.  66  ff ] ) . 

2.  An  Algebraic  Identity:  Given  the  discrete  sample  (at  equidistant  intervals) 

X1’  V  V  •  *  •’  V  •  •  •’  *T 

the  sample  covariance  at  lag  k  is  defined  as 


T-  |k| 

1 


(PT(k)  =  t  _  i 

0 


T  Li  xtxt+|k| 


for  k=  0,  +  1, .  .  .  ,  +  (T - 1)  * 
for  |k|>T 


v 


Recall  that  the  periodogram  is  given  by 


*T(a,)  =  2ffT 


x  e 
t 


icot 


t=  1 


*Note  that  for  large  lags  this  formula  for  cp( k)  yields  a  poor  estimate  of  y(k). 
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The  identity  in  question  reads 


*T^  2ttT 


T 

V 

>  x  e 

x- 

t=  1 


icat 


2tt 


|k|<T 


<^(k)e 


-  icak 


(1) 


regardless  of  assumptions  on  [12,  p.  52] .  Thus  4>rp(a>)  and  </>,p(k)  are  a  Fourier  pair. 
Now,  we  have  seen  that  for  a  zero  mean,  second  order  stationary  process  with  covariance 
sequence 


7(k)  =  <f 


Lxtxt+  Ik  | 


=  y(-k) 


a  spectral  density  exists  as 


k=-°° 


(2) 


Two  Fourier  pairs  —  one  for  the  sample,  one  for  the  process  —  are  now  in  evidence,  as 
shown  in  Table  1. 


Fourier  Transform  of  Function 

Sample 

$T<ca) 

Process 

f(co) 

y  (k) 

Table  1  Two  Fourier  Pairs 


Note  that  the  k)  provide  an  asymptotically  unbiased  estimate  of  the  y(k)  as  follows: 

T-  |k  |  T-  |k  | 

xtxt+|ki  MI 

t  =  1  t  =  1 


(f  [«^k)]  =(f 


<f 


XV  Ik  | 


T-  |k| 

=  ^  ^  y (k)  =(i.  -  )  y(k) 

t  =  1 


12 


p 

i.e. ,  (£^k)—  y(k);  the  sample  covariances  converge  in  probability*  to  the  process 
covariances.  Taking  expected  values  on  both  sides  of  (1),  we  have 


^  cp^. k)e 

|  k  |  <T 


=  k  I 

|  k  |  <T 

|k  |  <T 

Comparing  with  (2),  we  obtain 

Lim  <f  <i>  (w)  =  f  (on) 

T— 1 

Hence,  for  large  samples,  the  periodogram  appears  to  be  a  natural  estimator  of  the 
power  spectrum.  Surprisingly,  this  is  not  the  case.  It  has  been  found  by  experience 
that  a  harmonic  analysis  of  white  noise  (flat  spectrum)  produces  a  highly  spiked  spectrum 
(cf.  IV, B, 4).  The  reason  for  this  is  that  for  fixed  oo,  regardless  of  the  character  of  the 
noise, 

Lim  Var  4>^,(ajo)  =  2f  (a>Q), 

rp — ►GO 

whereas,  for  all  well  behaved  estimators  the  variance  should  go  to  zero  as  sample 
size  T— -°°.  It  can  in  fact  be  shown  that  for  Gaussian  processes  the  ratio  2$^(o;  )/f(o;  ) 
has  a  chi-square  distribution  (as  T-*-°°)  with  2  degrees  of  freedom.  This  result  stresses 
the  futility  of  appraising  the  significance  of  peaks  in  $^(0;)  even  as  T  increases,  when 

*  p 

t jT*  f  A  Lim  P{  |f  -f|>e}  =  0  ;  i.e.,  the  probability  that  differs  from  f  by  more 

T— 00 

than  e  goes  to  0  as  T  goes  to  00 . 
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the  spectrum  of  the  process  has  a  continuous  component  [2,  p.  4] . 

We  shall  turn  next  to  the  development  of  consistent  estimates  of  the 


spectrum . 

3.  Consistent  Estimates  (1):  Variance  Reduction:  As  an  estimate  of  f(u^), 

the  periodogram  has  a  non-zero  variance  independent  of  sample  size.  Bartlett  conceived 

the  idea  (1948)  of  using  a  variance  reducing  scheme  in  which  he  computes  his  estimate 

T 

of  the  power  spectrum  as  the  average  of  —  periodograms ,  each  computed  for  a  sub- 

m  2 

series  of  length  m.  If  the  variance  of  the  periodogram  estimate  of  f  (co  )  is  the 

m  9 

variance  of  the  new  estimate  may  be  expected  to  be  of  the  order  of  —  .  We  note 

here  that  as  more  segments  are  collected,  the  variance  is  reduced.  In  fact,  we  have 
a  consistent  estimate,  for  as  T-*°°  and  as  m  remains  constant,  var— '0.  This  fact  has 
been  the  starting  point  of  renewed  interest  in  the  theory  of  spectral  estimation.  We  shall 
see  later  that  a  price  is  paid  for  the  consistency  in  the  guise  of  an  asymptotic  bias. 

It  can  be  shown  that  Bartlett's  periodogram  averaging  scheme  is  equiv¬ 
alent  to  a  modification  of  the  classical  periodogram  computation  consisting  in  the  intro¬ 
duction  of  weighting  factors  for  the  sample  autocovariances.  Specifically,  mstead  of 
calculating 


27T 


I 

|  k  |  <T 


<p( k)e 

T 


-ito  k 
o 


as  an  estimate  of  f(wQ)>  one  uses 


where 


W = 


h  (k)  k)e 


|  k  |  <T 


icok 


r 

1- 


|k| 


h(k)=  V 


m 


0 


for  |  k  |  <.  m 
for  |k|>m 


In  this  manner,  an  entire  family  of  spectral  estimates  has  been  introduced  depending  on 
the  properties  of  the  weighting  function  h(k). 
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4.  Consistent  Estimates  (2):  Smoothing  Windows.  Although  is 

not  pointwise  consistent;  i.e.  ,  for  fixed  does  not  converge  in  probability*  to  f(cu)  as 
T— *>,  it  is  known  that  [12,  p.  58] 

X2 

P  f 

(o))do)—  \  f(u>)du>  , 

A  <  'i 


which  shows  that  the  averaged  spectrum  over  an  interval  can  be  consistently  estimated. 
The  left  side  may  be  modified  slightly  to  read 

OO 

f(w  )  =  \  W(o>)  , 

o  J  1 


with 


W(u)  =  S 


YX1 


for  X  <  u>  <  X_ 
1  2 


otherwise 


and  is  a  consistent  estimate  of  f(w  )  when  oj  is  in  the  interval  (X  ,  X  ).  It  is,  however, 

o  o  1  2 

biased,  since  it  estimates  an  averaged  or  smoothed  value  of  f(oj)  over  a  neighborhood 
of  u>o  on  the  u>  axis.  We  can  now  ask;  Is  there  a  way  in  which  we  can  tamper  with  W(- ) 
so  as  to  maintain  the  consistency  of  the  estimate  and  at  the  same  time  render  W(-) 
dependent  on  sample  size  T  in  such  a  way  as  to  reduce  the  bias  as  T-~°°  ? 

Heuristically,  we  started  with  W(-)  as  shown  in  Fig.  la, 


►  <*> 


U> 


Fig.  1.  Geometric  interpretation  of  the  periodogram  as  an  estimate  of 


*For  definition,  see  footnote  page  13. 
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and  estimated  the  average  of  f(co)  over  the  shaded  interval  in  Fig.  1(b).  Can  we  suitably 
modify  W(-)  so  that  it  depends  on  T  in  such  a  manner  as  to  become  increasingly  peaked 

as  T  increases  and  assume  the  shape  shown  in  Fig.  2,  all  the  while  satisfying  the  sub- 

00 

sidiary  condition  /  W  (co)dco  =  1?  If  this  program  can  be  carried  out,  the  values  of  f(u>) 

-00  1 

at  and  near  co^  will  be  given  more  and  more  weight  relative  to  values  of  f(<x>)  distant  from 


Fig.  2.  Smoothing  windows  sequence  and  their  dependence  on  sample  size. 

co  and  we  may  expect,  as  sample  size  increases,  to  reduce  the  bias  introduced  by 
averaging. 

5 .  Consistent  Estimates  (3):  Unified  Treatment.  Writing  the  periodogram 
4>^,(co)  as  the  Fourier  transform  of  the  sample  autocovariances  cp^ k)  we  have 

*T<U,)  =  U  I  vT(k)e 
|  k  |  <T 

In  Section  B.  3,  we  had  considered  estimates  tp(co)  of  the  form 

Vu>  =  -fa  Y,  hT<k>  «’T(k>e"lk“ 

|k|cT 

Let  the  Fourier  transform  of  h  (k)  be  represented  by 

V“>=s2  v*'*" 

|  k  |  <T 
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T  en,  the  convolution  theorem  applied  to  f_  gives 


7 r 

fT(^  =  2“  ^  V k)  y  k)e  =  j*  Ht  (w  -  X)  $T  (A)  dA  , 

J k  |  <T  -7 r 

which  demonstrates  the  equivalence  of  the  two  approaches  to  consistent  estimation 
<r  of  f(w).  ltp(k)  is  called  a  covariance  averaging  kernel  (lag  window),  is  called 

spectral  window.  Bartlett's  scheme  leads  directly  to  a  covariance  averaging  kernel. 

The  use  of  smoothed  spectra  leads  to  spectral  windows. 

6.  Statistical  Properties  of  Consistent  Estimates:  The  notion  of  bandwidth 
plays  an  important  role  in  the  formulation  of  the  statistical  properties  of  spectral  estimates 
of  the  type  we  have  just  considered.  The  bandwidth  of  a  spectral  window  W^,(-)  is 
defined  as  the  base  of  a  rectangle  which  has  the  same  area  and  the  same  maximum  height 
as  the  graph  of  W^,(-)  in  Fig.  2  (or  H_(*)  in  the  previous  section,  where  the  notation  H 
was  adopted  to  highlight  the  fact  that  h  and  H  are  a  Fourier  pair).  Formally, 

/  H„(o))dw 

-00  i 

BT  Max  |H  (w)  | 

Cl)  1 

where  the  dependence  on  sample  size  T  has  been  indicated.  We  would  like  B,_ -*■ 0  as 
T  —  °°,  to  improve  the  focusing  power  of  H^.( - )  on  the  peak  frequency  co^.  It  can  be 
shown  that  mean  square  error  and  bandwidth  are  related  as 

MSE  =  Bias2  +  Var  =  B  2p 

IT 

where  p  is  some  integer  >  1  (cf.  [10]  ).* 

The  important  point  is  to  note  the  mutually  antagonistic  effects  of  bias 

1 

and  variance.  Bandwidth  should  indeed  go  to  zero  as  T—  °°.  But  if  B  —  0  faster  than  —  , 
the  variance  term  will  become  very  large  indeed. 

*  p  is  the  order  of  the  highest  derivative  of  f(o;)  at  u;0.  Hence,  the  smoother  f(co)  at  c 
the  larger  is  p. 
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The  choice  of  a  specific  smoothing  window  or  kernel  averaging  sequence 
in  a  practical  situation  will  depend  on  how  sharp  an  estimate  of  f(k>Q)  one  desires;  with 
a  broad  bandwidth  the  variance  of  the  estimate  can  be  made  small;  on  the  other  hand,  if 
a  sharp  estimate  is  needed  (narrow  bandwidth,  small  bias)  the  variance  of  the  estimate 
increases.  We  might  further  note  that  if  the  spectrum  is  relatively  smooth  (changes 
slowly)  in  the  vicinity  of  then  its  average  value  over  an  interval  containing  w  should 
differ  but  little  from  f(u)o),  thus  insuring  a  small  bias  [20] . 
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III.  LIMITATIONS  INTRODUCED  BY  SAMPLING 


In  the  previous  chapter  we  have  concluded  that  a  periodic  regression  model  with 
known  base  frequency  is  of  limited  practical  value,  that  classical  periodogram  analysis 
leads  to  spurious  interpretations  if  a  continuous  component  is  present  in  the  spectrum 
and  that  covariance  averaging  kernels  (or  smoothing  windows)  are  of  decided  value  in 
estimating  spectra  for  processes  which  have  a  continuous  spectrum.  The  presentation 
was  made  in  terms  of  discrete  processes  because  this  is  usually  the  case  encountered 
in  practice;  but  the  statistical  properties  of  estimates  were  stated  in  terms  of  their 
asymptotic  behaviour.  Discrete  time-series  often  arise  from  sampling  continuous 
processes  and  even  in  the  case  of  discrete  processes  it  may  be  that  not  every  available 
piece  of  data  is  retained  in  the  analysis.  In  the  present  chapter,  we  will  briefly  dwell 
on  the  effects  of  sampling,  i.e. ,  of  retaining  only  a  finite  set  of  equi-distant  sample 
values.  Nor  will  the  asymptotic  character  of  the  results  of  the  previous  chapter  be 
strictly  exact  as  we  are  always  dealing  with  a  finite  sample;  they  serve  principally  as 
a  guide  in  the  choice  of  one  estimation  procedure  over  another. 

We  shall  be  dealing  here  with  a  specific  finite  sampled  sequence  of  values  rather 
than  with  an  ensemble  of  such  sequences.  We  shall  see  that  the  rejection  of  information 
inherent  in  sampling  introduces  ambiguities  and  limitations,  but  also,  that  in  a  profound 
sense  practically  all  the  information  in  the  sample  can  be  recovered. 

A.  What  is  Lost  in  Sampling  (1),  Uniqueness  of  Representation 

The  following  illustration  shows  clearly  the  ambiguities  which  are  introduced 
if  a  function  on  a  domain  is  sampled  at  a  subset  of  the  domain  of  definition.  In  brief, 
we  shall  "reconstruct"  a  function  from  its  sampled  values  in  two  different  ways  and  see 
what  happens  to  its  expansion.  Let  x(t)  be  defined  on  [0,T]  by 

x(t)  =  a  + 
o 

k=  1 

and  suppose  we  retain  only  values  of  x(t)  at  the  discrete  times  t  =  v ,  v-  1— T  . 


(a  .cos  kcot+  /3  sinkcot),  to  = 

K  Jv 


2tt 

T 
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3- PA-  2147 


(a)  ORIGINAL  FUNCTION  x(t) 


(b)  SAMPLED  VALUES  x 

v 


v  - 1/2  u  ♦  1/2 


(c)  STEP  FUNCTION  INTERPOLATION  x*(t)  (d)  POLYGONAL  INTERPOLATION  x  (t) 


Fig.  3*  Graphs  of  f(t)  and  its  approximations. 
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Then 


x  =  x(i')  ,  v  =  1  — T  ,  x  = 
v  o  T 

We  shall  reconstruct  the  Fourier  coefficients  of  x(t)  on  the  basis  of  the  sample  in  two 
different  manners: 

1.  By  assuming  that  x  is  reconstructed  as  a  step  function  (cf.  Fig.  3c) 

1  1 

x*(t)  =  x  for  v  <t<  y  +  — 
v  2  2 

2.  By  assuming  that  x  is  reconstructed  as  a  polygonal  line  connecting 
successive  ordinates 


x(t)  =  x  +  (x  -x  )  (t  -  v)  for  v<t<v+i  , 
v  v+ 1  v 


as  shown  in  Fig.  3d. 


x(t) 

x*(t) 

x(t) 

o 

< 

II 

o 

O 

< 

11 

*  o 
$ 

~a  =  A 
o  o 

< 

II 

w(V) 

K-\ 

vv  (0 

Table  2  Dependence  of  Fourier  Coefficients  on  Interpolating  Assumptions 
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^ 

The  least  squares  estimates  (a  ,  a,  ,  0, )  of  (a  ,  a,  ,  6, )  on  the  basis  of  the  sample 

okkokk 

are  given,  as  in  (II, A,  1)  by 

T 


A 

o 


X 

V 


V  =  1 


T 


x  coskwi' 

V 


2 

T 


x  sin  kcov 
v 

1 


A  straightforward  computation  yields  the  Fourier  coefficients  of  x*(t)  and  x(t)  in  terms 

of  the  (A  ,  A,  ,  B, ).  These  are  given  in  Table  2. 
o  k  k'  & 

Let  us  recapitulate.  We  have  reconstructed  a  function  from  its  sample  in  three 
different  ways:  (1)  as  a  trigonometric  series;  (2)  as  a  step  function;  (3)  as  a  polygonal 
line,  and  we  have  calculated  the  Fourier  coefficients  of  each  reconstruction.  It  is  now 
evident  that  these  depend  on  the  assumed  behaviour  of  the  function  between  the  sampled 
values . 


B.  What  is  Lost  in  Sampling  (2):  Folding. 

Suppose  a  set  of  discrete  observations  x^  are  obtained  as  a  result  of  sampling 
x(t)  at  fixed  intervals  AT.  The  frequency  components  of  x(t)  with  frequency  to  > 


AT 
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rad  /sec.  cannot  be  distinguished  from  those  with  frequency  to  in  the  range  (0,  — =-)  on 
the  sole  basis  of  the  sampled  values. 

To  see  this,  we  consider  first  the  case  of  a  single  sinusoid  and  show  that 


cos 


cos 


7 r 


VAT 


-o 


<P 
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for  all  t  =  k(AT),  i.  e. ,  for  all  integral  multiples  of  the  sampling  period  AT.  Indeed, 


kAT  4-  (p  =  k7r  + 


e  kAT+  (p 


(— 

Vat 


c  )  kAT  -  (p=  kir  - 


e  kAT+  c 


and  the  periodicity  and  symmetry  of  the  cosine  function  with  respect  to  k7r  complete 
the  proof. 

Figure  4  illustrates  the  symmetry  about  multiples  of  and  shows  why  this 
phenomenon  is  called  folding. 


Fig.  k.  Folding  illustrated. 


It  is  thus  seen  that  two  sinusoids  having  frequencies  symmetrically  placed  with  respect 
kir 

to  —  (k=  0,  +  1,  +  2,.  .  .)  cannot  be  distinguished  on  the  basis  of  values  sampled  at 
kAT.  This  has  been  called  "aliasing"  by  J.  W.  Tukey  in  the  sense  that  all  frequencies 
differing  from  a  given  frequency  u>  by  multiples  of^.,  {o>+  ^~}>are  indistinguishable 


AT 


from  go  the  so  called  principal  alias  of  cu,  which  is  that  frequency  among  all  {cu+ 

A  AT 

7 r 

which  falls  within  the  interval  0,  -^r 
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The  argument  just  developed  for  single  sinusoids  applies  equally  well  to  any 
function  x(t)  which  can  be  represented  as  a  superposition  of  sinusoids. 

In  general,  the  spectrum  f*(u;)  derived  from  a  sample  taken  at  rate  1/AT  is 
given  in  terms  of  the  true  spectrum  f(u>)  by 


OC 


k=  0 


7T 

The  folding  frequency  is  also  known  as  Nyquist  frequency.  A  graphical 

illustration  of  aliasing  is  shown  in  Fig.  5. 


3  -  PA  -  2 1 4  9 


Fig.  5-  Illustration  of  aliasing. 

7T  7T 

If  —  is  the  folding  frequency  of  the  sampling  process,  the  frequencies  in  f(ca)  beyond 

are  folded  back  and  added  to  the  original  graph  in  Fig.  5b.  A  spike  occurring  at  on ,  in  f(ca) 

S 

will  appear  in  f*(ca)  at  ca  .  ,  it's  principal  alias. 

S  A 

C.  What  is  Retained  in  Sampling:  The  Sampling  Theorem 

Admittedly,  sampling  discards  information.  There  is,  however,  one 
situation  in  which  nothing  is  lost. 
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Suppose  the  function  x(t)  contains  no  frequency  components  above 

W  cycles /second.  Then  x(t)  is  completely  determined  by  its  values  at  a  sequence 
.  .  1 

of  pomts  —  apart. 

In  particular,  if  the  sampled  values  are 
x(t)  can  be  uniquely  reconstructed  as 


(  2W  ) 


k  —  0,  +  1,  +  2, 


,  then 


+°° 


x(t)  = 


I 


k=  - 


sin  7r  (2Wt  -  k) 
7r  (2Wt  -  k) 


[19,33]  . 


Thus  a  band  limited  function  can  be  restored  completely  if  it  is  sampled  at  the  proper 

frequency.  Suppose  now  that  x(t)  is  not  band  limited,  i.  e. ,  there  are  frequency  com- 

1 

ponents  in  its  spectrum  beyond  W.  Then,  what  is  lost  by  sampling  at  a  rate  of  -r— 
is  precisely  what  has  been  described  as  aliasing,and  frequencies  beyond  W  are  folded 
back  into  the  interval  (0,  w^T).  If,  however,  x(t)  is  sampled  at  only  a  finite  number  of 
points,  say  for  k=0,  +  1,  +  2, .  .  . ,  +  n,  then  the  reconstruction  of  x(t)  as 

n 

x  ftv.V  sin ; (2Wt  -  k) 

iw  L  v2w  /  *  <2wt  -  k) 

k  =  -n 


is  but  one  of  infinitely  many  possible  reconstructions  agreeing  with  the  data  in  the 
sampled  interval.  Indeed,  for  any  sequence  y  ,  k  =  +  (n  +  1);  +  (n  +  2),  .  .  . 
the  function 


-(rH- 1)  <*> 
x2(t)  =  x1(t)  +  ^  +  / 

k=  -  oo  k=n+l 


sin  7 r  (2Wt  -  k) 
7r  (2Wt  -  k) 


agrees  with  x  (t)  for  k  =  0,  +  1,  +  2, .  .  . ,  +  n,  since  S^n  ^  vanishes  at  all 

1  —  _  _  7f  (2  wt  -  k) 

points  of  the  form  t=  ^ —  .  Thus,  values  sampled  at  equidistant  points  within  an  interval 

2W 

do  not  determine  values  outside  this  interval.  In  practice,  however,  [5]  empirical 
functions  are  always  band  limited,  nor  are  measurements  ever  so  precise  that  an  accurate 
reconstruction  could  be  obtained,  even  with  an  infinite  sequence  of  values.  In  this  sense, 
most  of  the  information  present  in  a  continuous  record  can  usually  be  extracted  from  a 
sample  taken  at  sufficiently  closely  spaced  ordinates. 
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D.  Implications  of  Sampling 

The  practical  implications  of  the  ambiguities  introduced  by  sampling  are 
several.  A  good  deal  depends  on  how  we  conceive  the  underlying  process  to  behave 
between  the  sampled  points.  If  it  is  nonexistent  there,  we  have  an  unlimited  choice; 
but  this  is  rarely  the  case.  In  general,  then,  aliasing  has  to  be  taken  into  account  in 
the  design  of  the  data  gathering  procedure  as  well  as  in  the  evaluation  of  the  results  of 
the  analysis. 

Since  the  sampled  spectrum  f*(w)  is  obtained  from  the  unsampled  spectrum 
f(w)  by  folding  all  frequencies  about  multiples  of  the  folding  frequency,  and  summing 
these  contributions  in  the  range  ^  0,  J,  we  must  insure,  if  we  wish  to  estimate  f(to) 
in  the  range  ^  0,-^p  j>  to  have  approximately  f(co)  =0  for  •  This  can,  in  general, 

be  done  by  filtering  the  data  through  a  low  pass  filter  before  sampling.  G.  M.  Jenkins  [14] 
has  suggested  that  in  any  practical  problem  there  are  3  frequencies  worth  considering: 

(1)  the  frequency  which  is  essentially  the  frequency  beyond  which 
the  recording  instruments'  frequency  response  function  is  negligible 
(e.g. ,  1-2%  of  its  maximum  value,  and  goes  to  0  asymptotically); 

(2)  the  frequency  beyond  which  there  is  barely  any  power  dissipated 
by  the  process  x(t)  ( ~  1  or  2%  of  the  total); 

(3)  a  good  a  priori  opinion  of  the  largest  frequency  of  interest  in  the 
study  of  the  process  at  hand. 

Frequency  has  a  transparent  significance,  as  it  limits  the  frequency  content  of  the 
data.  The  interplay  of  and  is  somewhat  more  subtle.  A  suitable  procedure  might 
be  to  base  the  sampling  rate  on  the  larger  of  and  If  >  u; well  and  good; 

0)3  =  it /AT,  so  that  AT  =  and  no  frequencies  beyond  those  of  interest  are  read; 

moreover,  there  is  negligible  aliasing  since  most  of  the  power  in  the  process  has  been 
dissipated  before  u ^  (<uO.  K  w-  >  w3’  comP°nents  of  f  beyond  could  be  folded 
relative  to  if  were  determining  the  sampling  rate.  Thus  we  select  to  determine 
AT  =  7r /^2’  providing  we  know  enough  about  the  process  to  choose  sensibly. 

One  of  the  costs  of  choosing  the  sampling  interval  on  the  basis  of  a  frequency 
higher  than  one  is  really  interested  in,  is  that  sampling  then  takes  place  at  a  much  higher 
rate  than  really  necessary,  giving  rise  in  turn  to  larger  quantities  of  data  to  process. 
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Another  cost  factor  is  related  to  the  nearness  of  the  frequency  of  interest  to 
the  folding  frequency;  the  closer  to— jr  >  the  harder  it  is  to  estimate,  in  the  sense  that 
a  larger  sample  is  required  (for  a  given  reliability  of  estimation)  than  for  a  frequency 
nearer  the  origin  of  the  frequency  axis.  Intuitively,  this  can  be  seen  in  the  following 
manner. 

Suppose  the  folding  frequency  is  u>^  (sampling  interval  AT  =— — —  )  and  a 

component  of  frequency  -e  (e  small)  is  present  in  the  signal.  The  sampling 

interval  corresponding  to  is  slightly  smaller  than  the  distance  between  successive 
zero  crossings  in  sin  cu^t, 


TT  PI  -  2I5Q~[ 


Fig.  6.  Ill  conditioning. 

and  ordinates  at  the  sampled  points  creep  up  extremely  slowly  to  the  peak  amplitude  as 
illustrated  in  Fig.  6.  Thus,  a  large  sample  may  be  needed  to  gain  confidence  in  the 
estimate  of  amplitude*.  As  a  rule  of  thumb,  it  is  usually  considered  unwise  to  estimate 
spectra  at  frequencies  exceeding  20  or  30%  of  the  folding  frequency,  as  the  reliability 
of  estimates  beyond  this  point  becomes  questionable.  There  is  a  great  deal  of  debate 
about  the  cutting  off  point,  depending  on  specific  situations  [  14,  28] . 


*The  nature  of  the  problem  is  one  of  ill  conditioning.  The  smaller  e,  the  more  difficult 
becomes  the  estimation.  An  analogy  might  clarify  matters  further:  it  is  easy  to  locate 
the  intercept  of  a  straight  line  if  its  slope  is  large.  It  becomes  much  more  difficult  to 
do  so  if  its  slope  is  very  small. 
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IV  MISCELLANEOUS  COMMENTS 


This  chapter  deals  with  several  points  which  could  not  properly  be  placed  in  the 
logical  sequence  of  the  previous  development  but  which  are  nevertheless  illuminating 
and  of  sufficient  practical  interest  to  be  included.  Some  of  these  are  direct  consequences 
of  the  convolution  theorem*  as  applied  in  (III ,  B,  5) ;  these  will  be  taken  up  first.  The 
remainder  of  the  chapter  will  present  some  useful  procedures  of  interest  in  actual 
computations,  and  some  examples  illustrating  the  methods  described  earlier. 

A.  The  Convolution  Theorem  and  Some  of  its  Consequences 


Let  us  recall  that  we  represented  the  periodogram 


(1) 


|k|<T 


as  the  Fourier  transform  of  the  sample  covariances,  and  that,  to  improve  its  usefulness 


as  an  estimator,  we  introduced  covariance  averaging  kernels  (lag  windows)  h^(k)  and 
obtained 


(2) 


|  k  |  <T 


as  an  estimate  of  f(u>).  Writing  the  Fourier  transform  of  h^(k)  as 


(3) 


|  k  |  <T 


we  noted  that,  by  the  convolution  theorem  applied  to  h^(k)  and 

7T 


(4) 


*The  convolution  theorem  states  that  if  the  following  are  Fourier  pairs  u(t)- — -  U(u>) 


OO 


-oo 
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We  are  now  in  a  position  to  clarify: 

-the  occurrence  of  negative  values  in  periodogram  estimates  of  the  spectrum 
although  the  latter  is  strictly  non -negative; 

-the  necessity  of  removing  D.  C.  components,  deterministic  components,  or 
time  dependent  non -zero  means  from  the  data  before  estimating  the  spectrum; 

-the  difference  in  appearance  between  plots  of  periodogram  and  smoothed  periodo¬ 
gram  respectively. 


1.  Negative  Values  of  the  Periodogram 

A  comparison  of  equations  (1)  and  (2)  shows  that  the  periodogram  (1)  can  be 
expressed  as  (2)  for  a  particular  choice  of  covariance  averaging  kernel,  namely 


I  1  for  |  k  |  <  T 
hT(k)  =  <J 

!  0  for  [k  |  >  T  . 

T 

Figure  7  illustrates  the  continuous  analog  of  (5)  and  its  Fourier  transform*  — 


(5) 


sin  ojT 
coT 


Fig.  7. 


A  ubiquitous  transform  pair. 


*There  are  some  analytical  complexities  introduced  by  the  "sampled  nature"  of  h^(k). 
We  do  not  go  into  full  detail,  but  the  argument  is  not  invalidated.  For  complete  details 
cf.  Blackmann  and  Tukey  [5] . 
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The  main  point  to  notice  is  that  sampling  perse  is  equivalent  to  the  introduction  of 
a  smoothing  window  H^cj)  which  assumes  negative  values,  so  that  the  spectrum 
estimate  (4)  allows  negative  weights  for  the  periodogram.  The  conclusion  is  that 

assumes  negative  values  with  non-zero  probability.  Large  negative  values  of 
the  periodogram  may  either  be  ignored  (e  g.  ,  if  it  is  felt  that  the  periodogram  is  well 
suited  to  the  modelling  assumptions)  or  they  may  be  taken  as  an  indication  that  other 
smoothing  and  estimation  schemes  should  be  attempted  in  order  to  ascertain  the 
stability  of  the  results. 

2.  The  Removal  of  D.C.  Components,  Time  Dependent  Non-Zero  Means, 

or  Deterministic  Components 

By  considering  the  effect  of  a  smoothing  window,  it  becomes  under¬ 
standable  that  D.C.  components  should  be  removed  before  processing.  Indeed,  every 
covariance  averaging  kernel  corresponds  to  a  window  with  some  equivalent  bandwidth 
over  which  the  said  D.C.  component  is  smeared  and  thus  made  to  contribute  to  the 
contamination  of  low  frequency  components.  A  similar  argument  may  be  made  for 
deterministic  components,  periodic  or  otherwise,  and  non -zero  means. 

Admittedly,  in  the  case  of  periodogram  calculations  the  bandwidth  of 
the  smoothing  window  can  be  made  arbitrarily  small  by  taking  a  large  enough  observa¬ 
tion  interval  (sample);  but  in  general,  since  it  is  desirable  to  use  covariance  averaging 
kernels  to  obtain  consistent  estimates,  it  is  advisable  to  remove  all  but  the  purely 
oscillatory  components  of  the  data. 

3.  Appearance  of  Periodogram  and  Smoothed  Periodogram  Plots 

The  shape  of  the  smoothing  window  corresponding  to  the  covariance 
averaging  kernel  used,  is  indicative  of  the  effect  of  smoothing  as  compared  to  the 
"unsmoothed"  estimate  of  the  spectrum.  If  a  spike  appears  in  the  unsmoothed  spectrum 
estimate,  the  smoothed  estimation  procedure  merely  redistributes  this  spike  over  the 
entire  axis  with  the  smoothing  window  as  envelope  for  the  spectral  components  at 
discrete  points.  Since  the  smoothing  window  is  generally  sharp  in  practice,  this 
redistribution  takes  place  mainly  over  a  small  neighborhood  of  the  frequency  at  which 
the  spectrum  is  estimated. 


31 


4. 


The  Effect  of  Non-Rational  Frequencies 
In  general,  when  the  sample  size  T  is  given,  spectral  analysis  con¬ 
sists  in  calculating  the  periodogram  (smoothed  or  unsmoothed)  at  a  discrete  set  of 

27T  i 

trial  frequencies.  u>  =  ,  j=  l  —  T/2,  under  the  assumption  that  x(t)  repeats  itself 

J  T 

exactly  after  T  observations,  apart  from  an  additive  error  term  (cf.  model  in  II,  A,  2). 
We  will  now  describe  the  effect  of  the  presence  of  a  frequency  component  whose  fre¬ 
quency  does  not  coincide  with  any  of  the  trial  frequencies. 

Let  a  process  be  represented  by  a  cyclical  deterministic  component, 
apart  from  additive  noise  ^ 


\  =  P  cos  ^  +  £t 


3-p«-?isT| 


2tt 
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where  k  is  the  multiple  of  —  nearest  to  —  X  (cf.  Fig.  8).  \ 


^(KH, 


2  IT 


(k+1) 


Fig.  8.  A  non -rational  frequency. 


It  can  be  shown  that  R,  .,  the  estimate  of  the  process  amplitude  p.  ,  (nonexistent!) 

k+j  r  k+j 

depends  asymptotically  (as  T  —  °°)  both  on  j  and  on  e  =  X-k  in  the  following  manner* 


k+j 


sin  7r  (j  -  e) 
Jr  (j  -c) 


so  that  the  irrationality  of  X  has  the  effect  of  smearing  the  correspondmg  frequency 
component  over  the  entire  spectrum.  The  presence  of  "irrational"  components  helps 
further  to  explain  why  periodogram  calculations  should  not  be  expected  to  yield  zero 


This  result  is  obtained  by  considering  the  Fourier  coefficients  of  cos  (o;  t+  ip)  on  a  base 
interval  which  is  not  an  integral  multiple  of  the  period.  In  the  worst  case,  when  e=  2 

€«  V4p2 


j  +  k 


~  2  (2j+l)2 

7T 
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estimates  for  amplitudes,  even  in  the  absence  of  periodic  components  at  the  trial 
frequency  (we  recall  that  another  source  of  non-zero  values  was  the  effect  of  the 
additive  random  noise  treated  in  the  model  of  Section  II,  A,  2). 

B.  Some  Useful  Procedures  and  Examples 

Several  procedural  details  are  offered  here  as  examples  of  the  kind  of 
improvements  which  can  be  brought  to  a  general  program  of  power  spectrum  estimation 
conceived  within  the  framework  of  the  theory  described  in  the  earlier  parts  of  this 
report.  Some  concrete  examples  are  provided  to  illustrate  some  of  the  developed  concepts. 
1.  Direct  Removal  of  Time  Dependent  Means 

In  the  processing  of  time  series  which  are  suspected  of  having  a 
non-zero  time  dependent  mean,  a  very  simple  scheme  exists  to  remove  this  component 
from  the  sample  covariances.  Instead  of  using  the  mean  lagged  product  as  an  estimate 
of  covariance,  calculate  instead 


(1) 


where  the  averages  are  taken  over  a  moving  strip  of  data. 

Actually,  this  procedure  assumes  that  the  process  x^  is  the  sum  of 
a  deterministic  component  m^  and  a  zero  mean  process  y  having  a  spectrum: 


Taking  expected  values  of  lagged  products,  we  obtain 


Now  formula  (1)  states  that  the  sample  autocovariance  <p(k)  of  the  zero  mean  process  y 
can  be  taken  as  the  average  lagged  product  of  the  sample  values  reduced  by  the  product 
of  the  sample  means.  Thus  the  spectral  analysis  proceeds  to  be  made  on  the  portion  of 
the  signal  for  which  the  theory  is  applicable. 

2.  Calculation  of  Sample  Covariances 

There  has  been  some  argument  that  the  inconsistency  of  the  periodogram 
as  an  estimate  of  the  power  spectrum  is  due  to  the  instability  of  the  sample  covariances 
as  estimates  of  the  process  covariances  as  the  lag  becomes  large.  In  practice,  there  are 
three  candidates  which  have  been  used  as  sample  covariances: 

T-  |k| 

=  kl  xtxt+iki'°ik-T-1 

t  =  1 

T-  |k| 

=  — -r— i  /  xx  |,|,0<k<T-l 

T-fk  |  Li  t  t+  |k| 

t  =  1 

T/2 

=  IS  xtxt+w-  °^r 

t=  1 

as  estimates  of  the  process  covariance  y(k).  We  can  compare  them  briefly  on  the  basis 
of  their  statistical  properties: 

k 

Bias:  cl  has  a  bias  equal  to  —  ,  which  goes  to  zero,  as 
T  —  °o,  for  every  fixed  k;  for  finite  samples  however,  the 
bias  persists.  c2  is  unbiased  for  any  lag  k  and  all  sample 
sizes  T. 

Positive  Definiteness:  cl  is  positive  definite,  which  is  a 
desirable  property  for  an  estimator  of  a  positive  quantity. 
c2  is  not  positive  definite. 


(cl)  (?(1)(k) 

T 

(c2)  (k) 

T 

(c3)  </>(3)(k) 

T 
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Mean  Square  Error:  As  an  estimate  of  yCk),  cl  has  a 

smaller  mean  square  error  than  c2  (cf.  [20]). 

As  l'ar  as  c3  is  concerned,  its  statistical  properties  have  not  been 
worked  out.  However,  as  an  estimate  of  y(k)  it  does  not  appear  quite  as  good  as  cl  and 
c2  since  it  uses  only  half  of  the  data  available.  But  this  in  itself  does  not  prove  that  it 
leads  to  poorer  estimates  of  the  power  spectral  density. 

3.  Examples 

The  importance  of  the  specific  choice  of  technique  is  illustrated  in 
Figs.  9,  10,  and  11.  Figure  9  is  a  plot  of  the  sample  autocovariance  of  a  set  of  data 
as  calculated  on  the  basis  of  formula  c3  of  Section  IV,  B,  1.  The  important  point  to  note 
is  that  two  periodic  components  are  clearly  present  in  the  autocovariance  (ACV)  function, 
one  of  relatively  high  frequency,  and  one  of  low  frequency. 

Figures  10  and  11  are  the  results  of  estimating  the  power  spectral 

density  in  two  different  ways.  The  method  used  to  obtain  Fig.  10  depends  on  formula  c2 

for  the  calculation  of  sample  autocovariances;  furthermore,  the  Fourier  transform  of 

111 

the  ACV  function  is  smoothed  with  Hanning  weights  ( — ,  — ,  — ). 

The  method  used  to  obtain  Fig.  11  proceeds  to  calculate  and  plot  the 
Fourier  transform  of  the  ACV  function  shown  in  Fig.  9. 

Note  that  to  a  hurried  observer,  the  crispness  of  Fig.  11  would 
recommend  itself  in  its  clear  identification  of  the  frequency  components  one  is  led  to 
expect  on  the  basis  of  a  glance  at  the  ACV  function  shown  in  Fig.  9.  A  more  seasoned 
analyst,  however,  will  recall  that  the  unsmoothed  periodogram  has  a  larger  mean  square 
error  than  the  smoothed  periodogram,  as  an  estimate  of  the  spectral  density. 

Admittedly,  some  crispness  could  be  lent  to  the  estimate  of  PSD  shown 
in  Fig.  10  by  omitting  the  final  smoothing  with  Hanning  weights.  Similarly,  smoothing 
with  Hanning  weights  would  reduce  the  sharpness  of  the  peaks  in  Fig.  11. 

In  conclusion,  it  must  be  said  that  there  is  at  present  no  single  optimal 
method  of  processing  data  when  little  is  known  apriori  about  either  the  process  generating 
the  data  or  the  characteristics  of  the  recording  device.  In  any  case,  it  is  desirable  to 
have  several  spectral  estimation  procedures  on  hand  to  gain  familiarity  with  new  data. 
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Autocovariance  f-unction  of  a  process 
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Fig.  10.  Estimate  of  power  spectrum  by  smoothed  (Hanning  weights)  periodogram. 
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Fig.  11.  Estimate  of  power  spectrum  after  removal  of  time  dependent  mean  (unsmoothed). 


4.  Illustration  of  Periodogram  Inconsistency  for  White  Noise 

The  probability  distributions  of  spectral  peaks  for  white  noise  are 
briefly  described  in  (II,  A,  2).  In  Section  (II,  B,  2)  it  was  further  mentioned  that  the 
periodogram  is  not  a  consistent  estimate  of  the  spectrum  of  a  zero  mean  second  order 
stationary  process,  and  that  as  a  consequence,  the  use  of  smoothed  spectra  came  to 
be  developed.  The  following  figures  (Figs .  12,  13,  and  14)  illustrate  the  situation 
rather  clearly. 

Figure  12  represents  a  sample  of  white  noise;  Fig.  13  is  a  plot  of 
its  autocorrelation;  and,  Fig.  14  is  a  plot  of  the  estimated  spectrum*.  The  main 
effect  is  to  observe  that  the  spectrum  is  highly  peaked  —  for  white  noise  —  even  in  the 
case  of  a  smoothed  periodogram,  as  shown  here.  For  a  classical  periodogram  of 
white  noise,  the  aspect  of  Fig.  14  would  be  more  jagged  still,  whereas  one  should  have 
expected  a  flat  spectrum. 


*The  spectrum  shown  here  is  a  smoothed  periodogram  with  a  Hanning  window.  A  classical 
periodogram  plot  was  not  available  at  the  time  of  writing. 
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Fig.  12.  A  sample  of  uncorrelated  white  noise. 
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Fig.  13-  Autocorrelation  function  of  white  noise. 
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Fig.  14.  Estimated  power  spectral  density  of  white 
noise  (smoothed  periodogram,  Hanning  weights). 
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V. 


BRIEF  OUTLINE  OF  A  PROGRAM  OF  TIME  SERIES  ANALYSIS 


The  basic  problem  of  time  series  analysis  may  be  formulated  as  a  question: 

‘'Given  an  empirical  time  series,  by  what  process  was  it 
generated? ,f  (P.  Whittle). 

This  is  a  familiar  question,  commonly  encountered  in  classical  statistical  estimation, 
and  usually  it  is  interpreted  as 

-  selecting  a  hypothetical  model  for  the  process 

-  setting  up  criteria  to  measure  the  degree  of  agreement  between  the  postulated 
model  and  the  observed  sequence. 

Given  two  models  to  choose  from,  we  can  always  decide,  according  to  our  established 
criteria,  which  is  the  more  adequate.  But  to  choose  a  reasonable  model  a  priori 
requires  some  wisdom.  In  the  present  context  (the  analysis  of  radar  data  from 
re-entering  vehicles),  it  would  appear  that  a  program  of  investigation  might  proceed 
as  described  below. 

First,  let  us  note  that  for  processes  with  a  discrete  spectrum  (e.  g.  ,  for  which 
the  cyclic  models  of  Section  (II,  A,  1  and  II,  A, 2)  are  applicable)  the  calculation  of 
smoothed  spectra  mars  only  slightly  the  occurrence  of  a  peak  by  redistributing  a 
portion  of  the  peak  amplitude  over  the  width  of  the  smoothing  window.  *  Hence,  in  the 
case  of  cyclical  processes,  little  is  lost  by  using  a  smoothed  periodogram  instead  of 
the  periodogram,  when  the  signal  to  noise  ratio  is  not  too  small.  On  the  other  hand, 
we  are  hardly  ever  faced  with  purely  cyclic  data  and  since  the  statistical  properties  of 
the  periodogram  are  so  poor  when  the  process  to  be  analyzed  is  oscillatory,  we  feel 
that  the  basic  spectral  estimation  tool  should  be  the  smoothed  periodogram.  ** 

*In  general,  only  equidistant  ordinates  are  considered  in  the  estimation  of  spectra.  If  a 
harmonic  term  falls  midway  between  two  ordinates,  its  amplitude  will  be  reduced  by  a 
4 

factor  of  roughly-*-  ^  •  41  at  these  ordinates,  in  the  case  of  classical  periodogram 

7T  L 

calculations.  This  may  or  may  not  be  enough  to  render  a  peak  non -significant  looking 
depending  on  the  signal  to  noise  ratio. 

**Since  the  estimation  will  take  place  via  covariance  averaging  kernels,  several  such 
kernels  corresponding  to  several  smoothing  windows  should  be  available  as  programmed 
options  in  any  computer  program  for  estimation  of  spectra. 

The  choice  of  a  smoothing  window  and/or  bandwidth  remains,  and  should  be  made  to 
depend  on  the  data  collection  scheme  and  the  resolution  and  precision  requirements 
pertinent  to  the  situation  under  investigation. 
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Second,  since  the  theory  is  applicable  only  to  zero  mean  second  order  stationary 
processes,  we  propose  to  use  the  techniques  of  spectral  estimation  only  on  the  residuals 
of  the  data,  after  removal  of  deterministic  components  and  non -zero  means. 

The  following  sequence  of  steps  is  envisaged. 

A.  The  Search  for  Regression 

As  a  first  step,  it  is  necessary  to  remove  deterministic  trends  from  the 
data.  Some  elegant  methods  exist  to  test  for  the  presence  of  polynomial  trends  adaptively 
on  the  degree  of  the  polynomial  (up  to  fixed  preassigned  maximum  degree)  against  the 
hypothesis  that  the  data  represents  pure  white  noise.  A  similar  adaptive  technique  can 
be  developed  for  other  families  of  regression  functions  if  it  should  appear  that  polynomials 
are  not  adequate  for  some  reason  or  other.  At  any  rate,  the  deterministic  component  is 
to  be  removed  from  the  data,  including  a  non-zero  mean  if  it  arises. 

B.  The  Search  for  Periodic  Components 

It  is  clear  that  a  periodogram  analysis  (smoothed  or  unsmoothed)  will  reveal 
the  presence  of  deterministic  cyclical  components  in  the  data.  Once  these  are  identified, 
they  can  be  removed  from  the  data.  For  identification  of  spurious  components,  see 
Section  C. 

C.  Spectral  Analysis 

At  this  point,  a  spectral  analysis  can  be  made,  using  a  variety  of  covariance 
averaging  kernels,  comparing  the  results  critically  among  themselves  and  particularly 
against  the  results  of  step  B  to  see  whether  the  removal  of  a  spurious  cyclic  component 
has  not  introduced  an  artificial  periodicity  in  the  data.  The  introduction  of  artificial 
periodicities  in  the  data  through  smoothing  is  far  from  being  an  academic  possibility, 
as  evidenced  by  a  theorem  due  to  Slutsky  (1927)  [11,  p.  95]  which  effectively  demon¬ 
strates  how  to  construct  a  linear  smoothing  scheme  which  introduces  an  arbitrarily 
selected  spurious  frequency. 

D.  Parameter  Estimation 

If  it  is  known  a  priori,  say  on  physical  grounds,  that  a  phenomenon  exhibits 
aspectrumof  a  given  functional  form  f(a>,  0 ^  .  .  . ,  0  )  depending  on  several  unknown 

parameters  0 0^,  .  .  .  ,  9^  these  parameters  may  be  estimated  on  the  basis  of  the 
empirical  time  series. 


*cf.  Reference  35. 
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It  should  be  possible  to  develop  a  testing  scheme  whereby  one  can  compare 
an  empirical  spectrum  to  an  a  priori  power  spectrum.  Such  a  technique  has  not  been 
found  in  the  literature  to  my  knowledge. 

E.  Data  Editing 

ri  here  is  a  preliminary  step  which  must  be  taken  before  any  estimation  and 
analysis  can  be  carried  out,  and  that  is  the  removal  of  spurious  data  and  of  debris  which 
usually  finds  its  way  past  any  recording  device.  Here,  however,  there  are  no  general 
prescriptions  as  every  case  must  be  treated  on  its  own  demerits.  For  an  example  of 
semi-automatic  data  editing  for  trajectory  data  the  reader  is  referred  to  Lincoln  Labora¬ 
tory  Report  PA-57  "A  Trajectory  Editing  and  Smoothing  Program  for  FPS-16  Radars” 
(LTnclassified)  where  some  of  the  many  possible  vicissitudes  of  data  editing  are  treated 
in  detail. 
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